Ce document présente un usage possible des packages R qui interfacent OpenStreetMap (OSM): géocoder une adresse, import de tuiles et de couches géographiques, import de points d’intérêt, routage et calcul d’itinéraires.
Il propose aussi des méthodes spatiales (carroyages, KDE, segments de lignes) et des visualisations associées pour restituer graphiquement la structure spatiale de semis de points ou de lignes.
En utilisant exclusivement des informations issues d’OSM, l’exercice consiste à identifier les secteurs accessibles à pieds dans le voisinage géographique d’un lieu défini par une adresse.
Reproduire l’analyse
Pour reproduire l’analyse, télécharger l’archive ZIP associée au dépôt https://github.com/riatecom/osm-elementr-2025-application, ouvrir le fichier proj.Rproj puis le fichier script.R qui reproduit en tout point le contenu du support de formation.
Cette analyse est reproductible en n’importe quel lieu, avec une adresse (adr) et un nom de ville (lib). Néanmoins les temps de calcul utilisant le serveur de démonstration d’OSRM ont été paramétrés pour ne pas le surcharger. L’exécution de ces blocs de code prennent du temps.
En conséquence, 6 jeux de données ont été préparés en amont à partir de l’adresse des gares SNCF de Lille, Toulouse, Grenoble, Nantes, Clermont-Ferrand et Compiègne, ainsi que de la station de métro Front Populaire. Les blocs de code relevant des calculs d’itinéraires / temps de trajets sont placés en commentaires dans le fichier script.R, mais sont fonctionnels.
lib adr
1 Grenoble 1 Place de la Gare, 38000 Grenoble, France
2 Clermont 40 Avenue de l'Union Soviétique, 63000 Clermont-Ferrand, France
3 Toulouse Gare SNCF Matabiau, 31500 Toulouse, France
4 Compiegne Gare de Compiègne, 60200 Compiègne, France
5 Condorcet Front Populaire, 93210 Saint-Denis, France
6 Nantes 27 Boulevard de Stalingrad, 44000 Nantes, France
7 Lille Lille Flandres, 59000 Lille, France
Nous présentons ici les résultats depuis la gare SNCF de Compiègne.
# 1. Packages qui interfacent des données OSM#install.packages('maposm', repos = 'https://riatelab.r-universe.dev')library(tidygeocoder) # géocodagelibrary(maptiles) # import de tuiles OSM (raster)library(osmdata) # import de données OSM (vecteur)library(maposm) # import de données OSM (couches géo)library(osrm) # calcul d'itinéraires# 2. Autres packages utilitaireslibrary(sf) # manipulation de données vectorielleslibrary(terra) # manipulation de données rasterlibrary(mapsf) # cartographie thématiquelibrary(maplegend) # légendeslibrary(spatstat) # analyse de semis de pointslibrary(mapview) # cartographie interactivelibrary(stplanr) # segmentiser des ligneslibrary(lwgeom)
1. Géocodage
Le package tidygeocoder(Cambon et al. 2025) permet d’utiliser des services de géocodage en ligne. Par défaut, le package utilise Nominatim qui se base sur les données OSM.
Passing 1 address to the Nominatim single address geocoder
Query completed in: 1 seconds
pt
Simple feature collection with 1 feature and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 2.82352 ymin: 49.42206 xmax: 2.82352 ymax: 49.42206
Geodetic CRS: WGS 84
# A tibble: 1 × 2
address geometry
* <chr> <POINT [°]>
1 Gare de Compiègne, 60200 Compiègne, France (2.82352 49.42206)
2. Import de tuiles OSM
Le package maptiles(Giraud 2025) permet de télécharger des tuiles chez différents fournisseurs, dont OpenStreetMap, au format raster.20 niveaux de zoom sont disponibles, selon la taille de l’objet géographique à cartographier. Un zoom de niveau 10 correspond à une échelle 1/500.000.
Plusieurs modèles de tuiles sont disponibles ici, certains nécessitent une inscription.
# Autre exemple de fournisseurtiles <-get_tiles(x = emprise,project =FALSE,crop =TRUE,provider ="OpenStreetMap.HOT",zoom =14)mf_raster(tiles, expandBB =c(rep(-.2,4)))mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_annotation(pt_reproj, txt =paste0("Gare de ", lib))mf_title("Tuiles OSM - OpenStreetmap.HOT")
3. Import de couches géographiques OSM
Pour la cartographie, le package maposm (pas sur le CRAN) se base sur les fonctions du package osmdata(Padgham et al. 2023) (partie suivante) pour télécharger les données OSM, au format vectoriel.
La fonction principale, om_get(), permet de récupérer certaines couches géographiques utiles (routes, parcs, bâti, etc) autour d’un point, dans un rayon donné. Le résultat prend la forme d’une liste d’objets sf.
res <-om_get(x =c(st_coordinates(pt)[1],st_coordinates(pt)[2]),r =2000)
mf_map(res$zone, col ="#f2efe9", border =NA, add =FALSE)mf_map(res$urban, col ="#e0dfdf", border ="#e0dfdf", lwd = .5, add =TRUE)mf_map(res$green, col ="#c8facc", border ="#c8facc", lwd = .5, add =TRUE)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5, add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$building, col ="#d9d0c9", border ="#c6bab1", lwd = .5, add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)mf_title("Autour de la gare")mf_arrow()mf_scale(size =500, scale_units ="m")mf_credits("OpenStreetMap contributors, 2025", pos ="bottomleft")
4. Import de tags avec osmdata
Le package osmdata permet d’extraire des données OSM à partir de l’API Overpass. Il nécessite dans un premier temps d’établir une bounding box au sein de laquelle seront extraites les données, puis de l’exécuter (fonction add_osm_feature()). Le résultat est une liste, qu’il faut ensuite convertir en objet sf.
Les objets OSM prennent la forme de nodes (points), ways (lignes ou polygones), ou relations (relation entre plusieurs objets). Ces éléments sont décrits par un ensemble de tags, recensées sous forme de clé-valeurs (key-values). Ici, nous recherchons dans la clé amenity les points (nodes) correspondant aux restaurants, bars et cafés (value = restaurant, bar, cafe). En dehors des clés-valeurs principales, d’autres attributs peuvent être présents, comme en l’occurence le type de cuisine, ou bien les horaires.
La fonction om_map() affiche les couches géographiques récupérées.
# Affichageom_map(res, title ="Autour de la gare")mf_map(poi_reproj, type ="typo", var ="amenity", cex = .7, pch =15, leg_pos ="topright", leg_val_cex = .9, leg_title_cex = .9, add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, lwd =2, col ="darkblue", add =TRUE)
Il est possible d’utiliser le package mapview(Appelhans et al. 2023) pour afficher les points ainsi que leurs attributs, de manière interactive.
mapview(res$zone, alpha.regions =0, layer.name ="Zone") +mapview(poi, layer.name ="Manger et boire")
osmextract
Le package osmextract(Gilardi et Lovelace 2024) permet d’extraire des données directement depuis une base de données OSM, sans avoir à passer par une API. Le plus petit niveau géographique auquel il est possible de télécharger les données correspond en France aux anciennes régions. Pour travailler sur Compiègne, il nous faudrait d’abord télécharger l’ensemble des données OSM de la Picardie, disponibles ici, et pesant plus de 240 Mb.
Le code ci-dessous effectue les mêmes manipulations que précédemment.
Une fois les points d’intérêt récupérés, il est possible de les agréger dans une maille.
# Création d'une grille régulière hexagonale de 300m de résolutiongrid <-st_make_grid(res$zone, cellsize =300, square =FALSE)grid <-st_intersection(grid, res$zone)grid <-st_sf(ID =1:length(grid), geom = grid)# Compter les POIs dans les carreauxinter <-st_intersects(grid, poi_reproj, sparse =TRUE)grid$n_poi <-lengths(inter)# Cartographieom_map(res)mf_map(grid, var ="n_poi", type ="choro", breaks =c(0,1,3,5,10,15),alpha = .6, border =NA, leg_title ="Nombre \nd'aménités",leg_val_rnd =0, leg_pos ="topright", add =TRUE)mf_title("Aménités gustatives - données carroyées")
6. Lissage (KDE)
Le package spatstat(Baddeley, Turner, et Rubak 2025) regroupe de nombreuses fonctions pour analyser des semis de points. Il génère des objets de classe ppp(planar point pattern), dans des fenêtres d’observation, de classe owin(observation window).
La fonction density.ppp() calcule des valeurs à partir d’une estimation par noyau (KDE). Elle prend en entrée une portée (sigma), ainsi qu’une résolution au sein de laquelle sera calculée la densité, en pixels (eps).
p <-as.ppp(X =st_coordinates(poi_reproj), W =as.owin(res$zone))ds <-density.ppp(x = p, sigma =100, eps =10, positive =TRUE)
Une fois les densités calculées, nous créons un objet de type SpatRaster. La densité étant calculée sur un mètre carré, nous la multiplions par 10.000 afin de récupérer une densité par hectare.
# Calcul densité de restaurants par hectarer <-rast(ds) *100*100crs(r) <-st_crs(poi_reproj)$wktr
class : SpatRaster
dimensions : 400, 400, 1 (nrow, ncol, nlyr)
resolution : 10, 10 (x, y)
extent : 312312.9, 316312.9, 6344781, 6348781 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857)
source(s) : memory
name : lyr.1
min value : 2.225074e-304
max value : 1.968529e+00
pal <-c( "#418e40", "#a5be6a", "#e2cf6a", "#c19849", "#ab6c32", "#842f1a")mf_map(res$zone, col ="#f2efe9", border =NA, expandBB =c(0,.3,0,0))mf_raster(r, leg_title ="N. restaurants / ha.", leg_pos ="left", pal = pal, add =TRUE)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5, add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(poi_reproj, pch =20, cex = .2, col ="black", add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)mf_title("La montagne culinaire - Compiègne")
Il est possible de convertir cette couche en format vecteur. Par défaut, la fonction as.polygons() arrondit les valeurs des pixels à l’entier.
# par défautr_poly <-as.polygons(r) |>st_as_sf()# tous les dixièmesr_poly <-as.polygons(r, round =TRUE, digits =1) |>st_as_sf()mf_map(res$zone, col ="#f2efe9", border =NA)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5,add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(r_poly, var ="lyr.1", breaks =quantile(r_poly$lyr.1),border =NA, lwd = .2, type ="choro", alpha =0.6,leg_title ="Densité\n (n/ha)", add =TRUE)mf_map(poi_reproj, pch =20, cex = .2, col ="black", add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)
Par aménité, discrétisation commune
7. Temps d’accès
Le package osrm(Giraud 2024) interface l’engin de routage OSRM et propose différentes fonctions pour le calcul de temps de trajet et d’itinéraires. Ces calculs s’effectuent sur la base de profils .lua qui pour chaque mode de transport définissent des algorithmes permettant la sélection de telle ou telle route, une vitesse par défaut selon le type de route, etc. Le profil voiture est disponible ici.
D’après la documentation de l’API OSRM, le nombre de requêtes sur le serveur de démo est limité à 512 par connection à l’API, espacées d’au moins 5 secondes. Pour des calculs plus volumineux (échelle nationale voire européenne), il est possible d’installer une instance OSRM sur son propre serveur (voir ici).
7.1. Depuis la gare
Pour calculer les temps de trajet (en minutes) et les distances (en mètres) entre la gare de Compiègne et l’ensemble des restaurants, bars et cafés situés dans un rayon de 2000m, nous utilisons la fonction osrmTable(), qui présente en sortie une liste composée de deux dataframes correspondant aux coordonnées des points d’origine et de destination, et d’un vecteur correspondant à la métrique calculée. Il est possible de calculer les distances et temps de trajet entre plusieurs couples origine / destination en une seule requête.
Nous construisons d’abord une matrice O/D composées de coordonnées X/Y des points d’origine et de destination.
Utilisation de l’API OSRM
Par souci de respect d’utilisation de l’API, une partie du code ci-dessous est passé en commentaire et non joué. Le résultat est enregistré sous format RDS pour être utilisé par la suite. Dans les morceaux de code suivants, le code apparaît dans des chunks sous option eval = FALSE.
# Création de la matrice O/Do <-data.frame(X =st_coordinates(pt)[, 1],Y =st_coordinates(pt)[, 2])row.names(o) <-1:nrow(o)d <-data.frame(X =st_coordinates(poi)[, 1],Y =st_coordinates(poi)[, 2])# 1 requête : 1 point d'origine, 92 points de destination#mat_dist <- osrmTable(src = o,# dst = d,# measure = "distance",# osrm.profile = 'foot')#saveRDS(mat_dist, "data/mat_dist.rds")# Intervalle de 5 secondes entre les requêtes#Sys.sleep(5)#mat_dur <- osrmTable(src = o,# dst = d,# measure = "duration",# osrm.profile = 'foot')#saveRDS(mat_dur, "data/mat_dur.rds")mat_dist <-readRDS("data/mat_dist.rds")mat_dur <-readRDS("data/mat_dur.rds")dist <-t(mat_dist$distances)dur <-t(mat_dur$durations)poi$dist <- distpoi$dur <- durapply(poi$dur, 2, summary)
1
Min. 1.000000
1st Qu. 7.600000
Median 10.500000
Mean 9.679747
3rd Qu. 11.750000
Max. 17.600000
apply(poi$dist, 2, summary)
1
Min. 73.0000
1st Qu. 570.5000
Median 789.0000
Mean 729.0759
3rd Qu. 881.5000
Max. 1501.0000
Il est possible d’atteindre 40 restaurants, cafés ou bars en moins de 11 min à pieds depuis la gare de Compiègne, et les trois quarts se situent à moins de 880m à pieds.
7.2. Depuis la grille
Afin de construire des indicateurs d’accessibilité, il est nécessaire de calculer les temps de trajet vers tous les restaurants / bars / cafés de la ville depuis chaque point de grille.
nrow(grid) *nrow(d)
[1] 14931
Nous avons en tout 14931 couples origine / destination. Pour éviter de surcharger le serveur de démo, nous découpons nos requêtes en autant de destinations. Ainsi, chaque requête calculera le temps de trajet entre nos 189 points d’origine et un point de destination.
Le code pour effectuer les calculs est présenté ci-dessous, mais a été joués en local. Les résultats sont enregistrés dans le dossier data.
osrmNearest()
Pour que l’engin de routage calcule un itinéraire, les points d’origine et de destination doivent être localisés sur une route. Si le point d’entrée ne se situe par sur une route, il sera déplacé sur le point le plus proche situé sur une route, en distance euclidienne. Cela peut avoir un impact sur l’itinéraire final ainsi que les temps de trajet, plus ou moins important selon le type d’espace (urbain ou rural).
# Coordonnées des points d'origine : centroïdes des cellules de la grilleo <-st_centroid(grid) |>st_transform(crs =4326)o$X <-st_coordinates(o)[, 1]o$Y <-st_coordinates(o)[, 2]o <-st_set_geometry(o[, c("X", "Y")], NULL)
Nous pouvons maintenant calculer plusieurs indicateurs d’accessibilité : temps de trajet au restaurant le plus proche, au deuxième, au troisième (question du choix), nombre de restaurants à moins de 5 minutes, 10 minutes, 15 minutes, etc.
mat <-read.csv("data/matfoot.csv")# Restaurant le plus proche pour chaque point de grillepoi_min <-aggregate(dur_foot ~ ID_ORIG, mat, FUN = min, na.rm =TRUE)# Récupérer l'identifiant de destination associépoi_min <-do.call(rbind,lapply(split(mat, mat$ID_ORIG),function(x) x[which.min(x$dur_foot), ]))# Cartographierpoi_min_sf <-merge(grid, poi_min,by.x ="ID", by.y ="ID_ORIG", all.x =TRUE)par(mfrow =c(1,1))om_map(res, title ="Temps d'accès au restaurant le plus proche")mf_map(poi_min_sf, var ="dur_foot", type ="choro", breaks =c(0,5,10,15,20,40),alpha = .6, border =NA, leg_title ="Temps\n(min.)",leg_val_rnd =0, pal ="Greens", leg_pos ="topright", add =TRUE)
Valhalla
Valhalla est un engin de routage open source basé sur les données OSM, au même titre que OSRM. Il fonctionne sur un système de tuiles contenant différents niveaux de routes, ainsi que les graphes routiers associés (les tuiles de niveau 0 contiennent les autoroutes, celles de niveau 2 les routes départementales). Les routes sont calculées au moment de la requête (tandis qu’elles sont pré-calculées sur OSRM) via différents algorithmes (A* pour les plus courts chemins, Contraction Hierarchies pour les longues distances). Valhalla propose de nombreux profils (vélo de location, multimodalité, taxi), ainsi que la possibilité de prendre en compte les pentes.
# moins de 5 minutesmat5 <- mat[mat$dur_foot <5 ,]mat5$n <-1mat5 <-aggregate(list(n_resto = mat5$n),list(ID_ORIG = mat5$ID_ORIG),FUN = sum)grid <-merge(grid[, "ID"], mat5,by.x ="ID", by.y ="ID_ORIG", all.x =TRUE)grid$n_resto[is.na(grid$n_resto)] <-0om_map(res, title ="Restaurants à moins de 5 minutes")mf_map(grid, var ="n_resto", type ="choro", breaks =c(0,1,5,10,30,max(grid$n_resto)),alpha = .6, border =NA, leg_title ="Nombre de\nrestaurants",leg_val_rnd =0, pal ="Reds", leg_pos ="topright", add =TRUE)
Le package osrm propose aussi des isochrones, à l’aide de la fonction osrmIsochrone().
isos <-st_read("data/mat.gpkg", layer ="isos", quiet =TRUE)isos <-st_transform(isos, crs ="EPSG:3857")isos <-st_intersection(isos, res$zone)om_map(res, title ="Isochrones autour de la gare de Compiègne")mf_map(isos, type ="typo", var ="isomax", alpha = .6, pal ="Geyser",border ="white", leg_pos =NA, add =TRUE)
8. Itinéraires
Pour aller plus loin, il est possible de récupérer les tracés des itinéraires, à l’aide de la fonction osrmRoute(). Attention, l’API route d’OSRM permet de calculer des itinéraires seulement entre deux points. Une requête correspond donc à un couple origine / destination.
Le package stplanr(Lovelace, Ellison, et Morgan 2024) permet initialement de modéliser et de travailler sur des systèmes de transport. Nous l’utilisons ici à des fins cartographiques, pour découper des lignes en segments.
routes <-st_read("data/mat.gpkg", layer ="routes", quiet =TRUE)routes <-st_transform(routes, crs ="EPSG:3857")routes$n <-1routes$distance <- routes$distance *1000# La fonction overline permet de compter le nombre d'occurence des tronçonsseg <-overline(routes, attrib ="n", fun = sum)routes <- routes[order(routes$distance, decreasing =TRUE), ]l_seg <-list()# La fonction line_segment découpe un objet sf LINESTRING en segments réguliersfor(i in1:nrow(routes)){ l_seg[[i]] <-line_segment( routes[i, ],segment_length =max(routes$distance) /20)}pal <-mf_get_pal(nrow(l_seg[[1]]), rev =TRUE, palette ="Blues 3")pal2 <-paste0(pal, "70")om_map(res, title ="Routes", theme ="dark")for(i innrow(routes):1){mf_map(l_seg[[i]], col = pal, lwd =1.6, add =TRUE)}leg(type ="cont", val =c(min(routes$distance), 1000, max(routes$distance)),pal = pal, pos ="left", title ="Distance (m)")
om_map(res, title ="Routes", theme ="dark")mf_map(seg, var ="n", type ="grad", col="#00f5f5", leg_pos =NA, breaks ="fisher", nbreaks =5, lwd =c(1,3,6,10,15), add =TRUE)mf_map(seg, var ="n", type ="grad", col="#fff70c", breaks ="fisher", nbreaks =5, lwd =c(.5,1.5,2.5,6,10), leg_val_rnd =0,leg_pos ="topright", leg_title ='N. Routes', add =TRUE)
om_map(res, title ="Routes", theme ="dark")mf_map(seg, var ="n", type ="grad", col="#00f5f570", leg_pos =NA, breaks ="fisher", nbreaks =5, lwd =c(2,5,9,14,18), add =TRUE)mf_map(seg, var ="n", type ="grad", col="#00f5f5", breaks ="fisher", nbreaks =5, lwd =c(1,3,7,11,13), leg_val_rnd =0,leg_pos ="topright", leg_title ='N. Routes', add =TRUE)
Bonus : la tournée des pizzeria
Je fais une escale d’une demie-heure à Compiègne et souhaite goûter toutes les pizzas de la ville en un minimum de temps depuis la gare. Avec la fonction osrmTrip(), c’est possible !
Cette fonction répond au problème du voyageur de commerce, qui consiste à déterminer le plus court circuit passant par chaque point une seule fois.
# Coordonnées des pizzeriapizz <- poi[!is.na(poi$cuisine) & poi$cuisine =="pizza" , "geometry"]# Ajouter le point de départ : la garepizz <-rbind(pt[1, "geometry"], pizz)pizz$id <-1:nrow(pizz)trip <-osrmTrip(loc = pizz,overview ="full",osrm.profile ="foot")trip <- trip[[1]]$tripst_write(trip, "data/mat.gpkg", layer ="trip", delete_layer =TRUE)
trip <-st_read("data/mat.gpkg", layer ="trip", quiet =TRUE)trips <-st_transform(trip, crs ="EPSG:3857")# Pour les labels : récupérer le début de la lignestart <- lwgeom::st_startpoint(trips)start <-do.call(rbind, start)start <-data.frame(id =rownames(trips), X = start[,1], Y = start[,2])start <-st_as_sf(start, coords =c("X", "Y"), crs ="EPSG:3857")par(mfrow =c(1,1))om_map(res, theme ="pizza", title ="La route de la pizza")mf_map(trips, col ="black", lwd =3, lty =3, add =TRUE)mf_map(start[start$id !=1, ], pch =20, cex =2, col ="black", add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_label(start, var ="id", cex =1.2, overlap =TRUE, pos =4, col ="black", halo =TRUE)
sum(trip$duration)
[1] 23.86833
D’après les données OSM, moins de 25 minutes suffisent pour faire la tournée des pizzas à Compiègne. On est large !!
Cambon, Jesse, Diego Hernangómez, Christopher Belanger, Daniel Possenriede, et Otto Hansen. 2025. tidygeocoder: An intuitive interface for getting data from geocoding services.https://CRAN.R-project.org/package=tidygeocoder.
Padgham, Mark, Bob Rudis, Robin Lovelace, Maëlle Salmon, et Joan Maspons. 2023. osmdata: Import ’OpenStreetMap’ Data as Simple Features or Spatial Objects. https://CRAN.R-project.org/package=osmdata.
Code source
---title: "Utiliser OpenStreetMap avec R"subtitle: "Mise en pratique"author: - name: Louis Laurian orcid: 0009-0002-7275-5536 affiliations: - ref: riate - name: Ronan Ysebaert orcid: 0000-0002-7344-5911 - name: Timothée Giraud orcid: 0000-0002-1932-3323affiliations: - id: riate name: Centre pour l’analyse spatiale et la géovisualisation (UAR RIATE) city: Paris state: France url: https://riate.cnrs.fr/date: 13 05 2025date-format: "D MMM YYYY"lang: frformat: html: embed-resources: true smooth-scroll: true fontsize: 0.9em code-tools: true code-fold: false toc: true toc-title: Utiliser OpenStreetMap avec R - Mise en pratique toc-depth: 3 toc-location: right bibliography: bib.bib linkcolor: "#d31329" css: custom.css include-after-body: "footer.html"editor_options: chunk_output_type: console---## ObjectifsCe document présente un usage possible des packages R qui interfacent OpenStreetMap (OSM): **géocoder une adresse**, **import de tuiles et de couches géographiques**, **import de points d'intérêt**, **routage** et **calcul d'itinéraires**.Il propose aussi des **méthodes spatiales** (carroyages, KDE, segments de lignes) et des **visualisations associées** pour restituer graphiquement la structure spatiale de semis de points ou de lignes. En **utilisant exclusivement des informations issues d'OSM**, l'exercice consiste à **identifier les secteurs accessibles à pieds dans le voisinage géographique d'un lieu défini par une adresse**. ## Reproduire l'analysePour reproduire l'analyse, télécharger l'archive ZIP associée au dépôt [https://github.com/riatecom/osm-elementr-2025-application](https://github.com/riatecom/osm-elementr-2025-application), ouvrir le fichier proj.Rproj puis le fichier script.R qui reproduit en tout point le contenu du support de formation. Cette analyse est reproductible en n'importe quel lieu, avec une adresse (adr) et un nom de ville (lib). Néanmoins les temps de calcul utilisant le serveur de démonstration d'OSRM ont été paramétrés pour ne pas le surcharger. L'exécution de ces blocs de code prennent du temps. En conséquence, 6 jeux de données ont été préparés en amont à partir de l'adresse des gares SNCF de Lille, Toulouse, Grenoble, Nantes, Clermont-Ferrand et Compiègne, ainsi que de la station de métro Front Populaire. Les blocs de code relevant des calculs d'itinéraires / temps de trajets sont placés en commentaires dans le fichier script.R, mais sont fonctionnels.```{r}cs <-read.csv("data/case_studies.csv", encoding ="UTF-8", sep =",")cs```Nous présentons ici les résultats depuis la **gare SNCF de Compiègne**. ```{r}cs <- cs[cs$lib =="Compiegne",]adr <- cs$adrlib <- cs$lib```## 0. Packages```{r, warning = FALSE, message = FALSE}# 1. Packages qui interfacent des données OSM#install.packages('maposm', repos = 'https://riatelab.r-universe.dev')library(tidygeocoder) # géocodagelibrary(maptiles) # import de tuiles OSM (raster)library(osmdata) # import de données OSM (vecteur)library(maposm) # import de données OSM (couches géo)library(osrm) # calcul d'itinéraires# 2. Autres packages utilitaireslibrary(sf) # manipulation de données vectorielleslibrary(terra) # manipulation de données rasterlibrary(mapsf) # cartographie thématiquelibrary(maplegend) # légendeslibrary(spatstat) # analyse de semis de pointslibrary(mapview) # cartographie interactivelibrary(stplanr) # segmentiser des ligneslibrary(lwgeom)```## 1. GéocodageLe package `tidygeocoder`[@R-tidygeocoder] permet d'utiliser des services de géocodage en ligne. Par défaut, le package utilise [Nominatim](https://nominatim.org/) qui se base sur les données OSM.```{r}pt <-geo(address = adr) |>st_as_sf(coords =c("long", "lat"), crs =4326)pt```## 2. Import de tuiles OSMLe package `maptiles`[@R-maptiles] permet de télécharger des tuiles chez différents fournisseurs, dont OpenStreetMap, au **format raster.** [20 niveaux de zoom](https://wiki.openstreetmap.org/wiki/Zoom_levels) sont disponibles, selon la taille de l'objet géographique à cartographier. Un zoom de niveau 10 correspond à une échelle 1/500.000.Plusieurs modèles de tuiles sont disponibles [ici](https://wiki.openstreetmap.org/wiki/Raster_tile_providers), certains nécessitent une inscription.```{r}emprise <-st_buffer(pt, 3000) |>st_bbox()tiles <-get_tiles(x = emprise,project =FALSE,crop =TRUE,cachedir ="cache",zoom =14)pt_reproj <-st_transform(pt, crs ="EPSG:3857")mf_raster(tiles, expandBB =c(rep(-.2,4)))mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_annotation(pt_reproj, txt =paste0("Gare de ", lib))mf_title("Tuiles OSM - Défaut")# Autre exemple de fournisseurtiles <-get_tiles(x = emprise,project =FALSE,crop =TRUE,provider ="OpenStreetMap.HOT",zoom =14)mf_raster(tiles, expandBB =c(rep(-.2,4)))mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_annotation(pt_reproj, txt =paste0("Gare de ", lib))mf_title("Tuiles OSM - OpenStreetmap.HOT")```## 3. Import de couches géographiques OSMPour la cartographie, le package `maposm` (pas sur le CRAN) se base sur les fonctions du package `osmdata`[@R-osmdata] (partie suivante) pour télécharger les données OSM, au **format vectoriel**.La fonction principale, `om_get()`, permet de récupérer certaines couches géographiques utiles (routes, parcs, bâti, etc) autour d'un point, dans un rayon donné. Le résultat prend la forme d'une liste d'objets sf.```{r}#| fig-width: 8#| fig-height: 8res <-om_get(x =c(st_coordinates(pt)[1],st_coordinates(pt)[2]),r =2000)names(res)mf_map(res$zone, col ="#f2efe9", border =NA, add =FALSE)mf_map(res$urban, col ="#e0dfdf", border ="#e0dfdf", lwd = .5, add =TRUE)mf_map(res$green, col ="#c8facc", border ="#c8facc", lwd = .5, add =TRUE)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5, add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$building, col ="#d9d0c9", border ="#c6bab1", lwd = .5, add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)mf_title("Autour de la gare")mf_arrow()mf_scale(size =500, scale_units ="m")mf_credits("OpenStreetMap contributors, 2025", pos ="bottomleft")```## 4. Import de tags avec osmdataLe package `osmdata` permet d'extraire des données OSM à partir de l'API [Overpass](https://overpass-api.de/). Il nécessite dans un premier temps d'établir une *bounding box* au sein de laquelle seront extraites les données, puis de l'exécuter (fonction `add_osm_feature()`). Le résultat est une liste, qu'il faut ensuite convertir en objet sf.Les objets OSM prennent la forme de *nodes* (points), *ways* (lignes ou polygones), ou *relations* (relation entre plusieurs objets). Ces éléments sont décrits par un ensemble de *tags*, recensées sous forme de clé-valeurs (*key-values*). Ici, nous recherchons dans la clé **amenity** les points (*nodes*) correspondant aux restaurants, bars et cafés (*value* = restaurant, bar, cafe). En dehors des clés-valeurs principales, d'autres attributs peuvent être présents, comme en l'occurence le type de cuisine, ou bien les horaires.```{r, warning = FALSE}# Définition d'une bounding box (emprise Compiègne)q <- opq(bbox = emprise, osm_types = "node")# Extraction des restaurants, bars et cafésreq <- add_osm_feature(opq = q, key = 'amenity', value = c("bar", "cafe", "restaurant"))poi <- osmdata_sf(req)poi <- poi$osm_points# Sélectionner les variablespoi <- poi[, c("osm_id", "name", "amenity", "cuisine")]# Intersectionpoi <- st_intersection(poi, st_transform(res$zone, crs = 4326))# Reprojectionpoi_reproj <- st_transform(poi , crs = "EPSG:3857")```La fonction `om_map()` affiche les couches géographiques récupérées.```{r}#| fig-width: 8#| fig-height: 8# Affichageom_map(res, title ="Autour de la gare")mf_map(poi_reproj, type ="typo", var ="amenity", cex = .7, pch =15, leg_pos ="topright", leg_val_cex = .9, leg_title_cex = .9, add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, lwd =2, col ="darkblue", add =TRUE)```Il est possible d'utiliser le package `mapview`[@R-mapview] pour afficher les points ainsi que leurs attributs, de manière interactive.```{r}mapview(res$zone, alpha.regions =0, layer.name ="Zone") +mapview(poi, layer.name ="Manger et boire")```::: callout-note### osmextractLe package `osmextract`[@R-osmextract] permet d’extraire des données directement depuis une base de données OSM, sans avoir à passer par une API. Le plus petit niveau géographique auquel il est possible de télécharger les données correspond en France aux anciennes régions. Pour travailler sur Compiègne, il nous faudrait d'abord télécharger l'ensemble des données OSM de la Picardie, disponibles [ici](https://download.geofabrik.de/europe/france/picardie.html), et pesant plus de 240 Mb.Le code ci-dessous effectue les mêmes manipulations que précédemment.```{r, eval = FALSE}library(osmextract)osm_pt <- oe_get(place = "Picardie", layer = "points", extra_tags = "amenity", quiet = TRUE)osm_pt <- osm_pt[osm_pt$amenity %in% c("bar", "cafe", "restaurant"), ]osm_pt <- st_transform(osm_pt, crs = "EPSG:3857")osm_pt <- st_filter(osm_pt, res$zone, .predicates = st_intersects)```:::## 5. CarroyageUne fois les points d'intérêt récupérés, il est possible de les agréger dans une maille.```{r}#| fig-width: 8#| fig-height: 8# Création d'une grille régulière hexagonale de 300m de résolutiongrid <-st_make_grid(res$zone, cellsize =300, square =FALSE)grid <-st_intersection(grid, res$zone)grid <-st_sf(ID =1:length(grid), geom = grid)# Compter les POIs dans les carreauxinter <-st_intersects(grid, poi_reproj, sparse =TRUE)grid$n_poi <-lengths(inter)# Cartographieom_map(res)mf_map(grid, var ="n_poi", type ="choro", breaks =c(0,1,3,5,10,15),alpha = .6, border =NA, leg_title ="Nombre \nd'aménités",leg_val_rnd =0, leg_pos ="topright", add =TRUE)mf_title("Aménités gustatives - données carroyées")```## 6. Lissage (KDE)Le package `spatstat`[@R-spatstat] regroupe de nombreuses fonctions pour analyser des semis de points. Il génère des objets de classe **ppp** *(planar point pattern)*, dans des fenêtres d'observation, de classe **owin** *(observation window)*.La fonction `density.ppp()` calcule des valeurs à partir d'une estimation par noyau (KDE). Elle prend en entrée une portée (*sigma*), ainsi qu'une résolution au sein de laquelle sera calculée la densité, en pixels (*eps*).```{r}#| fig-cap: "@moraga2023"#| echo: falseknitr::include_graphics(rep('img/kde.png'))``````{r}p <-as.ppp(X =st_coordinates(poi_reproj), W =as.owin(res$zone))ds <-density.ppp(x = p, sigma =100, eps =10, positive =TRUE)```Une fois les densités calculées, nous créons un objet de type *SpatRaster*. La densité étant calculée sur un mètre carré, nous la multiplions par 10.000 afin de récupérer une densité par hectare.```{r, warning = FALSE}#| fig-width: 8#| fig-height: 7# Calcul densité de restaurants par hectarer <- rast(ds) * 100 * 100crs(r) <- st_crs(poi_reproj)$wktrpal <- c( "#418e40", "#a5be6a", "#e2cf6a", "#c19849", "#ab6c32", "#842f1a")mf_map(res$zone, col = "#f2efe9", border = NA, expandBB = c(0,.3,0,0))mf_raster(r, leg_title = "N. restaurants / ha.", leg_pos = "left", pal = pal, add = TRUE)mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)mf_map(poi_reproj, pch = 20, cex = .2, col = "black", add = TRUE)mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)mf_title("La montagne culinaire - Compiègne")```Il est possible de convertir cette couche en format vecteur. Par défaut, la fonction `as.polygons()` arrondit les valeurs des pixels à l'entier.```{r}#| fig-width: 8#| fig-height: 8# par défautr_poly <-as.polygons(r) |>st_as_sf()# tous les dixièmesr_poly <-as.polygons(r, round =TRUE, digits =1) |>st_as_sf()mf_map(res$zone, col ="#f2efe9", border =NA)mf_map(res$water, col ="#aad3df", border ="#aad3df", lwd = .5,add =TRUE)mf_map(res$railway, col ="grey50", lty =2, lwd = .2, add =TRUE)mf_map(res$road, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(res$street, col ="white", border ="white", lwd = .5, add =TRUE)mf_map(r_poly, var ="lyr.1", breaks =quantile(r_poly$lyr.1),border =NA, lwd = .2, type ="choro", alpha =0.6,leg_title ="Densité\n (n/ha)", add =TRUE)mf_map(poi_reproj, pch =20, cex = .2, col ="black", add =TRUE)mf_map(res$zone, col =NA, border ="#c6bab1", lwd =4, add =TRUE)```Par aménité, discrétisation commune```{r, echo = FALSE}# Cafés et bars ensemblepoi_reproj$amenity[poi_reproj$amenity %in% c("bar", "cafe")] <- "bar-cafe"p <- list()ds <- list()r <- list()sel <- levels(as.factor(poi_reproj$amenity))for(i in 1:length(sel)){ p[[i]] <- as.ppp( X = st_coordinates(poi_reproj[poi_reproj$amenity == sel[i] ,]), W = as.owin(res$zone)) ds[[i]] <- density.ppp(x = p[[i]], sigma = 100, eps = 10, positive = TRUE) r[[i]] <- rast(ds[[i]]) * 100 * 100 crs(r[[i]]) <- st_crs(poi_reproj)$wkt r[[i]] <- as.polygons(r[[i]], round = TRUE, digits = 1) |> st_as_sf()}# Récupérer l'amplitudeamp <- do.call(rbind, r)amp <- sort(c(unique(amp$lyr.1)))bks <- quantile(amp)pal <- mf_get_pal(n = length(bks) -1, palette = "Reds", alpha = .6, rev = TRUE)par(mfrow = c(1,2))mf_map(res$zone, col = "#f2efe9", border = NA)mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)mf_map(r[[1]], var = "lyr.1", breaks = bks, type = "choro", border = NA, leg_title = "Densité\n (n/ha)", leg_pos = "bottomleft", pal = pal, add = TRUE)mf_map(poi_reproj[poi_reproj$amenity == "bar-cafe" ,], pch = 20, cex = .2, col = "black", add = TRUE)mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)mf_title("Densité de bars et cafés")mf_map(res$zone, col = "#f2efe9", border = NA)mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)mf_map(r[[2]], var = "lyr.1", type = "choro", border = NA, breaks = bks, leg_title = "Densité\n (n/ha)", pal = pal, leg_pos = "bottomleft", add = TRUE)mf_map(poi_reproj[poi_reproj$amenity == "restaurant" ,], pch = 20, cex = .2, col = "black", add = TRUE)mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)mf_title("Densité de restaurants")```## 7. Temps d'accèsLe package `osrm`[@R-osrm] interface l'engin de routage [OSRM](https://project-osrm.org/) et propose différentes fonctions pour le calcul de temps de trajet et d'itinéraires. Ces calculs s'effectuent sur la base de **profils .lua** qui pour chaque mode de transport définissent des algorithmes permettant la sélection de telle ou telle route, une vitesse par défaut selon le type de route, etc. Le profil voiture est disponible [ici](https://github.com/Project-OSRM/osrm-backend/blob/master/profiles/car.lua).D'après la [documentation de l'API OSRM](https://project-osrm.org/docs/v5.24.0/api/#), le nombre de requêtes **sur le serveur de démo** est limité à **512** par connection à l'API, espacées d'au moins **5 secondes**.Pour des calculs plus volumineux (échelle nationale voire européenne), il est possible d'installer une instance OSRM sur son propre serveur (voir [ici](https://rcarto.github.io/posts/build_osrm_server/)).### 7.1. Depuis la garePour calculer les temps de trajet (en minutes) et les distances (en mètres) entre la gare de Compiègne et l'ensemble des restaurants, bars et cafés situés dans un rayon de 2000m, nous utilisons la fonction `osrmTable()`, qui présente en sortie une liste composée de deux dataframes correspondant aux coordonnées des points d'origine et de destination, et d'un vecteur correspondant à la métrique calculée. **Il est possible de calculer les distances et temps de trajet entre plusieurs couples origine / destination en une seule requête.**Nous construisons d'abord une matrice O/D composées de coordonnées X/Y des points d'origine et de destination.::: callout-warning#### Utilisation de l'API OSRMPar souci de respect d'utilisation de l'API, une partie du code ci-dessous est passé en commentaire et non joué. Le résultat est enregistré sous format RDS pour être utilisé par la suite. Dans les morceaux de code suivants, le code apparaît dans des chunks sous option *eval = FALSE.*:::```{r}# Création de la matrice O/Do <-data.frame(X =st_coordinates(pt)[, 1],Y =st_coordinates(pt)[, 2])row.names(o) <-1:nrow(o)d <-data.frame(X =st_coordinates(poi)[, 1],Y =st_coordinates(poi)[, 2])# 1 requête : 1 point d'origine, 92 points de destination#mat_dist <- osrmTable(src = o,# dst = d,# measure = "distance",# osrm.profile = 'foot')#saveRDS(mat_dist, "data/mat_dist.rds")# Intervalle de 5 secondes entre les requêtes#Sys.sleep(5)#mat_dur <- osrmTable(src = o,# dst = d,# measure = "duration",# osrm.profile = 'foot')#saveRDS(mat_dur, "data/mat_dur.rds")mat_dist <-readRDS("data/mat_dist.rds")mat_dur <-readRDS("data/mat_dur.rds")dist <-t(mat_dist$distances)dur <-t(mat_dur$durations)poi$dist <- distpoi$dur <- durapply(poi$dur, 2, summary)apply(poi$dist, 2, summary)```Il est possible d'atteindre `r round(nrow(d)/2, -1)` restaurants, cafés ou bars en moins de `r ceiling(median(poi$dur))` min à pieds depuis la gare de Compiègne, et les trois quarts se situent à moins de `r round(apply(poi$dist, 2, summary)[5],-1)`m à pieds.### 7.2. Depuis la grilleAfin de construire des indicateurs d'accessibilité, il est nécessaire de calculer les temps de trajet vers tous les restaurants / bars / cafés de la ville depuis chaque point de grille. ```{r}nrow(grid) *nrow(d)```Nous avons en tout `r nrow(grid) * nrow(d)` couples origine / destination. Pour éviter de surcharger le serveur de démo, nous découpons nos requêtes en autant de destinations. Ainsi, chaque requête calculera le temps de trajet entre nos 189 points d'origine et un point de destination. Le code pour effectuer les calculs est présenté ci-dessous, mais a été joués en local. Les résultats sont enregistrés dans le dossier *data*.::: callout-important#### osrmNearest()Pour que l'engin de routage calcule un itinéraire, les points d'origine et de destination doivent être localisés sur une route. Si le point d'entrée ne se situe par sur une route, il sera déplacé sur le point le plus proche situé sur une route, en distance euclidienne. Cela peut avoir un impact sur l'itinéraire final ainsi que les temps de trajet, plus ou moins important selon le type d'espace (urbain ou rural).:::```{r, warning = FALSE}# Coordonnées des points d'origine : centroïdes des cellules de la grilleo <- st_centroid(grid) |> st_transform(crs = 4326)o$X <- st_coordinates(o)[, 1]o$Y <- st_coordinates(o)[, 2]o <- st_set_geometry(o[, c("X", "Y")], NULL)``````{r, eval = FALSE}# Les points de destination restent les mêmesmatfoot <- list()for(i in 1:nrow(d)){ matfoot[[i]] <- osrmTable(src = o[, c("X", "Y")], dst = d[i , c("X", "Y")], measure = "duration", osrm.profile = 'foot') Sys.sleep(5) matfoot[[i]]$ID_ORIG <- c(1:nrow(o)) matfoot[[i]]$ID_DEST <- rep(i, nrow(o)) matfoot[[i]] <- data.frame("ID_ORIG" = matfoot[[i]]$ID_ORIG, "ID_DEST" = matfoot[[i]]$ID_DEST, "duration" = matfoot[[i]]$durations) colnames(matfoot[[i]]) <- c("ID_ORIG", "ID_DEST", "dur_foot")}matfoot <- do.call(rbind, matfoot)write.csv(matfoot, "data/matfoot.csv", row.names = FALSE)```Nous pouvons maintenant calculer plusieurs indicateurs d'accessibilité : temps de trajet au restaurant le plus proche, au deuxième, au troisième (question du choix), nombre de restaurants à moins de 5 minutes, 10 minutes, 15 minutes, etc.```{r}#| fig-width: 8#| fig-height: 8mat <-read.csv("data/matfoot.csv")# Restaurant le plus proche pour chaque point de grillepoi_min <-aggregate(dur_foot ~ ID_ORIG, mat, FUN = min, na.rm =TRUE)# Récupérer l'identifiant de destination associépoi_min <-do.call(rbind,lapply(split(mat, mat$ID_ORIG),function(x) x[which.min(x$dur_foot), ]))# Cartographierpoi_min_sf <-merge(grid, poi_min,by.x ="ID", by.y ="ID_ORIG", all.x =TRUE)par(mfrow =c(1,1))om_map(res, title ="Temps d'accès au restaurant le plus proche")mf_map(poi_min_sf, var ="dur_foot", type ="choro", breaks =c(0,5,10,15,20,40),alpha = .6, border =NA, leg_title ="Temps\n(min.)",leg_val_rnd =0, pal ="Greens", leg_pos ="topright", add =TRUE)```::: callout-important#### ValhallaValhalla est un engin de routage open source basé sur les données OSM, au même titre que OSRM. Il fonctionne sur un système de **tuiles** contenant différents niveaux de routes, ainsi que les graphes routiers associés (les tuiles de niveau 0 contiennent les autoroutes, celles de niveau 2 les routes départementales). Les routes sont calculées au moment de la requête (tandis qu'elles sont pré-calculées sur OSRM) via différents algorithmes (A* pour les plus courts chemins, *Contraction Hierarchies* pour les longues distances). Valhalla propose de nombreux profils (vélo de location, multimodalité, taxi), ainsi que la possibilité de prendre en compte les pentes.:::```{r}#| fig-width: 8#| fig-height: 8# moins de 5 minutesmat5 <- mat[mat$dur_foot <5 ,]mat5$n <-1mat5 <-aggregate(list(n_resto = mat5$n),list(ID_ORIG = mat5$ID_ORIG),FUN = sum)grid <-merge(grid[, "ID"], mat5,by.x ="ID", by.y ="ID_ORIG", all.x =TRUE)grid$n_resto[is.na(grid$n_resto)] <-0om_map(res, title ="Restaurants à moins de 5 minutes")mf_map(grid, var ="n_resto", type ="choro", breaks =c(0,1,5,10,30,max(grid$n_resto)),alpha = .6, border =NA, leg_title ="Nombre de\nrestaurants",leg_val_rnd =0, pal ="Reds", leg_pos ="topright", add =TRUE)```Le package osrm propose aussi des isochrones, à l'aide de la fonction `osrmIsochrone()`.```{r, eval = FALSE}isos <- osrmIsochrone(loc = pt, breaks = seq(0,45,5), res = 60, osrm.profile = "foot")st_write(isos, "data/mat.gpkg", layer = "isos", delete_layer = TRUE)``````{r, warning = FALSE}#| fig-width: 8#| fig-height: 8isos <- st_read("data/mat.gpkg", layer = "isos", quiet = TRUE)isos <- st_transform(isos, crs = "EPSG:3857")isos <- st_intersection(isos, res$zone)om_map(res, title = "Isochrones autour de la gare de Compiègne")mf_map(isos, type = "typo", var = "isomax", alpha = .6, pal = "Geyser", border = "white", leg_pos = NA, add = TRUE)```## 8. ItinérairesPour aller plus loin, il est possible de récupérer les tracés des itinéraires, à l'aide de la fonction `osrmRoute()`. Attention, l'API route d'OSRM permet de calculer des itinéraires seulement entre deux points. **Une requête correspond donc à un couple origine / destination.**```{r, eval = FALSE}routes <- list()o <- data.frame(X = st_coordinates(pt)[, 1], Y = st_coordinates(pt)[, 2])row.names(o) <- 1for(i in 1:nrow(d)){ routes[[i]] <- osrmRoute(src = o, dst = d[i,], overview = "full", osrm.profile = 'foot') Sys.sleep(5)}routes <- do.call(rbind, routes)st_write(routes, "data/mat.gpkg", layer = "routes", delete_layer = TRUE)```Le package `stplanr`[@R-stplanr] permet initialement de modéliser et de travailler sur des systèmes de transport. Nous l'utilisons ici à des fins cartographiques, pour découper des lignes en segments.```{r, warning = FALSE}routes <- st_read("data/mat.gpkg", layer = "routes", quiet = TRUE)routes <- st_transform(routes, crs = "EPSG:3857")routes$n <- 1routes$distance <- routes$distance * 1000# La fonction overline permet de compter le nombre d'occurence des tronçonsseg <- overline(routes, attrib = "n", fun = sum)routes <- routes[order(routes$distance, decreasing = TRUE), ]l_seg <- list()# La fonction line_segment découpe un objet sf LINESTRING en segments réguliersfor(i in 1:nrow(routes)){ l_seg[[i]] <- line_segment( routes[i, ], segment_length = max(routes$distance) / 20)}pal <- mf_get_pal(nrow(l_seg[[1]]), rev = TRUE, palette = "Blues 3")pal2 <- paste0(pal, "70")om_map(res, title = "Routes", theme = "dark")for(i in nrow(routes):1){ mf_map(l_seg[[i]], col = pal, lwd = 1.6, add = TRUE)}leg(type = "cont", val = c(min(routes$distance), 1000, max(routes$distance)), pal = pal, pos = "left", title = "Distance (m)")om_map(res, title = "Routes", theme = "dark")mf_map(seg, var = "n", type = "grad", col= "#00f5f5", leg_pos = NA, breaks = "fisher", nbreaks = 5, lwd = c(1,3,6,10,15), add = TRUE)mf_map(seg, var = "n", type = "grad", col= "#fff70c", breaks = "fisher", nbreaks = 5, lwd = c(.5,1.5,2.5,6,10), leg_val_rnd = 0, leg_pos = "topright", leg_title = 'N. Routes', add = TRUE)om_map(res, title = "Routes", theme = "dark")mf_map(seg, var = "n", type = "grad", col= "#00f5f570", leg_pos = NA, breaks = "fisher", nbreaks = 5, lwd = c(2,5,9,14,18), add = TRUE)mf_map(seg, var = "n", type = "grad", col= "#00f5f5", breaks = "fisher", nbreaks = 5, lwd = c(1,3,7,11,13), leg_val_rnd = 0, leg_pos = "topright", leg_title = 'N. Routes', add = TRUE)```## Bonus : la tournée des pizzeriaJe fais une escale d'une demie-heure à Compiègne et souhaite goûter toutes les pizzas de la ville en un minimum de temps depuis la gare. Avec la fonction `osrmTrip()`, c'est possible !Cette fonction répond au [problème du voyageur de commerce](https://fr.wikipedia.org/wiki/Probl%C3%A8me_du_voyageur_de_commerce), qui consiste à déterminer le plus court circuit passant par chaque point une seule fois.```{r, eval = FALSE}# Coordonnées des pizzeriapizz <- poi[!is.na(poi$cuisine) & poi$cuisine == "pizza" , "geometry"]# Ajouter le point de départ : la garepizz <- rbind(pt[1, "geometry"], pizz)pizz$id <- 1:nrow(pizz)trip <- osrmTrip(loc = pizz, overview = "full", osrm.profile = "foot")trip <- trip[[1]]$tripst_write(trip, "data/mat.gpkg", layer = "trip", delete_layer = TRUE)``````{r}#| fig-width: 8#| fig-height: 8trip <-st_read("data/mat.gpkg", layer ="trip", quiet =TRUE)trips <-st_transform(trip, crs ="EPSG:3857")# Pour les labels : récupérer le début de la lignestart <- lwgeom::st_startpoint(trips)start <-do.call(rbind, start)start <-data.frame(id =rownames(trips), X = start[,1], Y = start[,2])start <-st_as_sf(start, coords =c("X", "Y"), crs ="EPSG:3857")par(mfrow =c(1,1))om_map(res, theme ="pizza", title ="La route de la pizza")mf_map(trips, col ="black", lwd =3, lty =3, add =TRUE)mf_map(start[start$id !=1, ], pch =20, cex =2, col ="black", add =TRUE)mf_map(pt_reproj, pch =24, cex =1.3, col ="darkblue", lwd =2, add =TRUE)mf_label(start, var ="id", cex =1.2, overlap =TRUE, pos =4, col ="black", halo =TRUE)sum(trip$duration)```D'après les données OSM, moins de 25 minutes suffisent pour faire la tournée des pizzas à Compiègne. On est large !!## [](img/logo.png)